Biological materials.

Adult Lamellibrachia lumysi specimens were collected from seep localities in the Mississippi Canyon at 754 m depth in the Gulf of Mexico (N 28°11.58’, W 89°47.94’), using the R/V Seward Johnson and Johnson Sea Link in October 2009. All samples were frozen at 80˚C following collection.

Genome sequencing and assembly.

Vestimentum tissue was dissected from one individual of the worm, and high molecular weight genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer’s protocols. Sequencing a total of six paired-end or mate-pair genomic DNA libraries with insert sizes ranging from 180 bp to 7 kb were performed by The Genomic Services Lab at the Hudson Alpha Institute in Huntsville, Alabama on an Illumina HiSeq 2000 platform (see details in Table S1). Paired-end libraries (180 bp, 400 bp, 750 bp) were prepared using the 125 bp TrueSeq protocols, and mate-pair libraries (3-5 kbp, 5-7 kbp) were generated using the Illumina Nextera Mate Pair Library Kit followed by size selection. Moreover, a 10X sequencing library was constructed using the 10X Chromium protocol (10X genomics) at the Hudson Alpha Institute. The finished library was sequenced on an Illumina HiSeqX platform, using paired 151 bp reads with a single 8 bp index read.

Our workflow of genome assembly was shown in Fig. S1. The paired-end and 10X raw reads were checked using FastQC v0.11.5 (Andrews and others 2010) and quality filtered (Q score >30) using Trimmomatic v0.36 (Bolger, Lohse, and Usadel 2014). The estimation of genome size, level of heterozygosity and repeat contents of the Lamellibrachia genome was determined by analysing the kmer histograms generated from the paired-end libraries using Jellyfish v2.2.3 (Marçais and Kingsford 2011) and GenomeScope (Vurture et al. 2017) (Fig. S2). The Mate-pair reads were trimmed and sorted using NxTrim v0.3.1 (O’Connell et al. 2015) which can recognize and trim the artificial Nextera mate-pair circulation adapters. Only reads from category “mp” (true mate-pair reads) and “unkonwn” (mostly large insert size reads) were used for downstream scaffolding anlaysis. Reads from “pe” (paired-end reads) and “se” (single ends) categories were discarded.

Given that high heterozygosity of Lamellibrachia genome, all reads were assembled using Platanus v1.2.4 (Kajitani et al. 2014) with a kmer size of 32. Scaffolding was conducted by mapping Illumina paired-end and mate-pair reads to contigs generated by Platanus using SSPACE v3.0 (Boetzer and Pirovano 2014). Gaps in the scaffolds were then filled with GapCloser v1.12 (Luo et al. 2012). Redundant allele scaffolds were further removed using Redundans v0.13c with default settings (Pryszcz and Gabaldón 2016). Genome assembly quality was assessed using QUAST v3.1 (Gurevich et al. 2013). Completeness of obtained genome was assessed using BUSCO v3(Waterhouse et al. 2017) with Metazoa_odb9 database (978 busco genes).

Transcriptome assembly and analysis

Total RNA was extracted from the plume, vestimentum and trophosome tissue from the same individual of Lamellibrachia specimen using Trizol. RNA-Seq of adult tissues from plume, vestimentum, and trophosome was performed using Illumina HiSeq 2000 platform in Hudson Alpha. After quality checking and trimming of raw sequencing reads, transcripts were assembled de novo with Trinity. Transcript isoforms with high similarity (≥ 95%) were removed with CD-HIT-EST v4.7 (Li and Godzik 2006). Transcript abundance was estimated with Bowtie v2.2.9 (Langmead and Salzberg 2012) and RSEM v1.2.26 (Li and Dewey 2011) by mapping reads back to the transcriptomic assembly based on transcripts per million (TPM). A tissue-specifically expressed gene was defined as a gene that had over 75% of the total transcripts in a particular tissue based on TPM (Albertin et al. 2015). GO enrichment analysis was performed based on the GO annotation using PANNZER2 (Törönen, Medlar, and Holm 2018). Statistically overrepresented GO terms of trophosome-specific genes were identified using Lamellibrachia gene models as background with AgriGO (Du et al. 2010).

Proteomics characterization.

Proteomics analyses of Lamellibrachia trophosome tissue were performed by Proteomics & Metabolomics facility at Colorado State University. Briefly, samples were transferred to 5mL tubes appropriate for use in the Bullet Blender Storm 5 (Next Advanced). Freshly cleaned, 3.2mm stainless steel beads were added at an approximate 1:1 sample: bead volume ratio (1-5beads). The tissue/bead mixture was then combined with 1X PBS (Hyclone) spiked with 1X Halt Protease Inhibitor cocktail (Thermo Scientific) for an approximate 2:1 liquid: solid ratio (150-600μl). Tissue homogenization was achieved using speed setting 12 for 5 minutes. The stainless steel beads were removed and 0.1mm glass beads added at approximately 1:1 sample:bead ratio. Further homogenization was achieved at speed 8 for 3 minutes. Homogenate was then transferred to 1.5mL microcentrifuge tubes and subjected to cup horn sonication (Amplitude 70, 10s pulse, 20s off, 7 minutes total exposure time; Qsonica) while suspended in ice water. Insoluble material was pelleted by brief centrifugation at 3000xg. The resulting supernatant was subjected to protein quantification by the Pierce BCA Protein Assay Kit (ThermoFisher-Pierce) following manufacturer’s instructions. Absorbance of reacted samples was measured at 550nm and protein quantification stemmed from a quadratic fit of the Bovine Serum Albumin (BSA) standard curve.

50 μg total protein was aliquoted from each sample and processed for in-solution trypsin digestion as previously described (Schauer et al. 2013). A total of 0.5μg of peptides were then purified and concentrated using an on-line enrichment column (Waters Symmetry Trap C18 100Å, 5μm, 180 μm ID x 20mm column). Subsequent chromatographic separation was performed on a reverse phase nanospray column (Waters, Peptide BEH C18; 1.7μm, 75 μm ID x 150mm column, 45°C) using a 90 minute gradient: 5%-30% buffer B over 85 minutes followed by 30%-45%B over 5 minutes (0.1% formic acid in ACN) at a flow rate of 350 nanoliters/min. Peptides were eluted directly into the mass spectrometer (Orbitrap Velos, Thermo Scientific) equipped with a Nanospray Flex ion source (Thermo Scientific) and spectra were collected over a m/z range of 400–2000, positive mode ionization, using a dynamic exclusion limit of 2 MS/MS spectra of a given m/z value for 30 s (exclusion duration of 90 s). The instrument was operated in FT mode for MS detection (resolution of 60,000) and ion trap mode for MS/MS detection with a normalized collision energy set to 35%. Compound lists of the resulting spectra were generated using Xcalibur 3.0 software (Thermo Scientific) with a S/N threshold of 1.5 and 1 scan/group.

Tandem mass spectra were extracted, charge state deconvoluted and deisotoped by ProteoWizard MsConvert (version 3.0). Spectra from all samples were searched using Mascot (Matrix Science, London, UK; version 2.6.0) against gene models of Lamellibrachia host and symbiont genomes (derived from (Y. Li, Liles, and Halanych 2018)) assuming the digestion enzyme trypsin. Mascot was searched with a fragment ion mass tolerance of 0.80 Da and a parent ion tolerance of 20 PPM. Oxidation of methionine and carbamidomethyl of cysteine were specified in Mascot as variable modifications. Search results from all samples were imported and combined using the probabilistic protein identification algorithms (Keller et al. 2002) implemented in the Scaffold software (version Scaffold_4.8.4, Proteome Software Inc., Portland, OR) (Searle, Turner, and Nesvizhskii 2008). Peptide thresholds were set such that a low peptide FDR was achieved based on hits to the reverse database (Käll et al. 2007) or Protein Prophet (Nesvizhskii et al. 2003). Protein identifications were accepted if they could be established at greater than 99.0% probability and contained at least 1 identified peptide. Protein probabilities were assigned by the Protein Prophet algorithm (Nesvizhskii et al. 2003). Proteins that contained similar peptides and could not be differentiated based on MS/MS analysis alone were grouped to satisfy the principles of parsimony.

Genome annotation.

Our genome annotation workflow was shown in Fig. S3. Gene models of Lamellibrachia genome were constructed following the Funannotate pipeline 1.3.0 (https://github.com/nextgenusfs/funannotate). Briefly, repeptitive regions in the Lamellibrachia genome were identified using RepeatModeler v1.0.8 (Smit and Hubley 2008) and were subsequently soft-masked using RepeatMasker v4.0.6 (Chen 2004). RNA-Seq data from different tissue were leveraged to improve the accuracy of gene prediction. RNA-Seq data were assembled de novo into transcriptomes using Trinity v2.4.0 (Haas et al. 2013) and HISAT 2.1.0 (Kim, Langmead, and Salzberg 2015) was used to algin RNA-Seq reads to the Lamellibrachia assembly. Tanscriptome assemblies were then passed to PASA pipeline v2.3.3 (Haas et al. 2003) to identify high-quality gene models. The aligned RNA-Seq data was used to train the ab initio gene predictions using AUGUSTUS v3.3 (Stanke et al. 2006). Protein alignments from the SwissProt database to Lamellibrachia assembly was generated using exonerate (Slater and Birney 2005) and Trinity/PASA transcripts were aligned to the genome using Minimap2 v2.1 (H. Li 2018). The tRNA genes were identified using tRNAscan-SE v1.3.1 (Lowe and Eddy 1997). Finally, EvidenceModeler 1.1.0 (Haas et al. 2008) was used to combine all the evidences of gene prediction from protein alignments, transcript alignments, and ab initio predictions to construct high-quality gene models. Finally, functional annotations of predicted gene models were analyzed using several curated databases. KEGG orthology was assigned using the KEGG Automatic Annotation server. Gene models were further annotated with domain structure and protein identity by InterProScan (Zdobnov and Apweiler 2001) and SwissProt database, respectively. Secreted proteins were predicted using SignalP (Petersen et al. 2011) and Phobius (Käll, Krogh, and Sonnhammer 2007) using InterProScan.

Gene family analysis

After all-to-all Diamond v1.0 (Buchfink, Xie, and Huson 2014) BLASTP searches against 22 selected lophotrochozoan proteomes (Table S3), orthology groups (OGs) were identified using Orthofinder with a default inflation parameter (I=1.5). Gene ontology annotation was performed using PANTHER v13.1 (Mi et al. 2016) with the PANTHER HMM scoring tool (pantherScore2.pl). Gene family expansion and contraction was estimated using CAFÉ v2.1 (De Bie et al. 2006). For each gene family, CAFÉ generated a family-wide P value, with a significant P value indicating a possible gene-family expansion or contraction event. A branch-specific P value was also generated for each branch/node using the Viterbi method. A family-wide P value less than 0.01 and a branch/node Viterbi P value less than 0.001 was considered as a signature of gene family expansion/contraction for a specific gene family and specific species, respectively, as suggested in manual.

In addition to the above gene family expansion or contraction analysis, the protein-coding sequences were used for discovery of genes under positive selection in from L. lumysi in comparison with selected lophotrochozoan species (including Capitella teleta, Lottia gigantea, Helobdella robusta, Octopus bimaculoides, Crassostrea virginica, Notospermus geniculatus and Phoronis australis) using PosiGene, which is an automated pipeline for genome-wide detection of positively selected genes , with default settings (p < 0.05). In this positive selection analysis through the PosiGene pipeline, L. lumysi was used as the reference and anchored species for ortholog assignment of each OG.

Manual annotation of gene families with potential interest.

In addition to the automatic annotation process mentioned above, we manually annotated genes that are general interest, such as hemoglobin gene families, genes related to amino aicd synthesize, or immunology. These genes were specifically selected a priori based on our experience after a revision of the available publications in the field. For hemoglobin genes, datasets used for HB gene family analyses, identified homologs, and their respective accession numbers are available in Fig. S# .Hbs and linker sequences of interest were obtained from “Lamellibrahcia” genome and assembled siboglinid transcriptomes (see molecular clock analysis above) via diamond blastp (evalue cutoff 1E-5) by utilizing Hb and linker sequences acquired from GenBank of siboglinids as queries (Figure #). Sequences with best hists to target proteins were annotated for protein domain architecture usig the Pfam databases included in InterProscan. After manual removal of redundant and incorrect sequences (e.g., sequences are too short or lack of globin domain), we used MAFFT 7.2.15 (Katoh and Standley 2013) to align HB amino acid sequences. Maximum likelihood analyses were performed in IQTree under the best-fitting models for associated partition schemes determined by Modelfinder implemented in IQTree with ultrafast bootstrapping of 1000 replicates.

For genes related to immunity,MIKE, I

Genes related to amino acids synthesis were from Lamellibrachia, Lamellibrachia symbionts, C. teleta were performed by using the KEGG2 KAAS genome annotation web server and then visualized by the KEGG Mapper Reconstruct Pathway too. A BLASTp search of protein sequences from genome annotation against the Swiss-Prot database was used to search for proteins that were missing in the visualized KEGG pathway.

Phylogenomics and Molecular clock analysis

Analysis of siboglinid phylogeny was conducted from publically available 16 transcriptomic data in conjunction with Lamellibrachia proteome (Table S2). Sequence assembly, annotation, homology evaluation, gene tree construction, parsing of genes trees to OGs, and supermatrix construction were conducted with Agalma (Dunn, Howison, and Zapata 2013). The final analyses presented here contained 13 siboglinids and four outgroups based on current understanding of annelid phylogeny. Maximum likelihood analyses were performed in IQTree under the best-fitting models for associated partition schemes determined by Modelfinder implemented in IQTree 1.6.3 (Nguyen et al. 2014) with ultrafast bootstrapping of 1000 replicates.

For molecular clock analysis, a relaxed molecular clock with a lognormal distribution and a Yule tree model was used in BEAST 2 v2.5.1 (Bouckaert et al. 2014). A calibration was placed on the node representing the most recent common ancestor (MRCA) of Osedax using a normal distribution with a mean of 100 MYA and a standard deviation of 10 following the findings of (Danise and Higgs 2015). Another calibration was placed on the node of MRCA of Serpulida and Sabellida using a normal distribution with a mean of 267 MYA (Sanfilippo et al. 2017). Molecular clock analyses with BEAST 2 consisted of two independent runs with 1 million MCMC generations sampled every 1000 generations. Convergence was checked and confirmed by comparing trace plots in Tracer making sure the effective sample size of each parameter was greater than 100 and that stationarity appeared to have been achieved. Log and tree files were combined using Logcombiner. A maximum clade credibility tree with mean heights was calculated using TreeAnnotater. The resulted time-calibrated tree was plotted using R package, phyloch, strap, and OutbreakTools. Bayesian inference using a molecular clock resulted in identical branching patterns as analysis with IQTree. ## Results and Discussion

Genome features

The genome of Lamellibrachia was sequenced using a combination of Illumina paired-end, mate-pair and 10x techniques. Results from high-throughput sequencing and genome assembly for Lamellibrachia are presented in Table S3 with at least 300 fold coverage using a combination of Illumina paired-end, mate-pair and 10X genomic sequencing with a variety of insert sizes. The haploid genome assembly size is ~687 Mb, similar to the result of genome size estimation using GenomeScope (~646 Mb). The N50 value of the assembled scaffolds and contigs is 373 Kb and 24 Kb, respectively. Although N50 lengths and assembly quality of Lamellibrachia are comparable to those of other annelids (e.g. Capitella teleta, Hebdella robusta), the overall genome completeness measured by BUSCO (~ 95%) is one of the highest among other lophotrochozoans (Table S4), indicating the completeness of the genome assembly. With the support of RNA-seq data from different tissue, we estimated Lamellibrachia genome contains 38,998 gene models. The genome also exhibits a similar level of heterozygosity (0.6%) and repetitive contens (36.92%), compared to other lophotrochozoan genomes. The repeat contents of Lamellibrachia genome was shown in Table S#.

Nutritional adaptations to a symbiotic life stage

Figure 2. Lamellibrachia lacks many genes related to amino acid biosynthesis. (A) Genes present in Lamellibrachia symbiont and Capitella teleta genomes but absent in Lamellibrachia genome. (B) Example of Lysine biosynthesis pathway in Lamellibrachia. Red boxes represent missing genes in the pathway.

Adult Lamellirbachia worms lack a digestive tract and rely solely on endosymbionts for nutrition. Noticeably, only 57 genes associated with amino acid biosynthesis were found in the Lamellibrachia host genome, of which eight were identified as proteins (Fig S#). In contrast, the genome of the Lamellibrachia symbiont (110 genes) contains essentially complete gene sets for the biosynthesis of all 20 proteinogenic amino acids and of 11 vitamins and cofactors. The lack of 30 genes in “Lamellibrachia” compared to “Capitella” genome (90 genes) are involved in the biosynthesis of 13 out of 20 amino acids (e.g. five key enzymes in the Lysine biosynthesis pathway Fig. #). As amino acids are indispensable for protein biosynthesis in the host, the sparse presence of amino acid synthesis-related genes in Lamellibrachia genome may indicate that the host depends on its symbionts for supply with amino acids and cofactors. Moreover, we found a large gene expansion of nutrient uptake ABC transports in Lamellibrachia compared with other lophotrochozoans (Table #). These findings are consistent with previous biochemical analyses suggest that Riftia is also dependent on its bacterial symbiont for the biosynthesis of polyamines that are important for its metabolism and physiology (Minic and Hervé 2003). Intrestingly, a similar mechanism of host defficiency of amino acid synthesis has also been suggested in thiotrophic vesicomyid symbioses of Calyptogena magnifia (Newton, Girguis, and Cavanaugh 2008) and vent musssel symbioses of Bathymodiolus azoricus (Ponnudurai et al. 2017). These results emphasize the host deficiency of amino acid syntheses and exchanges of metabolites between two partners might to be a common strategy for chemosynthetic symbioses.

In addition to the immediate release of fixed carbon and provide amino acids by the symbionts, there is evidence of a second possible nutritional mode - digestion of symbionts directly by the host, as showed by the detection of abundant host-derived digestive enzymes in the trophosome tissue (Table S#). Although previous observations indicated that symbiont could be digested by Riftia (Bright, Keckeis, and Fisher 2000), the direct evidence and genetic mechanism related to symbiont digestion have not been characterized. We identified fifteen the host proteins related to lysosomal proteases that were highly expressed and detected as proteins in the trophosome tissue, such as Saposin and multi copies of Cathepsin (Table S#). Lysosomes contain an array of digestive enzymes used for protein degradation in eukaryotic cells and are thought to play an essentail role in symbiont digestion of B. azoricus (Ponnudurai et al. 2017). We furthermore identified 19 major components of proteasome as proteins in the trophoseome tissue, indicating a potential role in protein degradation of symbiont digestion. Thus, our results suggest that lysosomal host proteases and proteasome may facilitate the degradation of symbiont proteins during symbiont digestion, and maintaining population size of symbionts within trophosome, or possibly both.

Because the trophosome tissue of Lamellibrachia contains sulfur-oxidizing symbionts, we also identified ~ 200 bacterial proteins in the trophosome tissue to gain insights into the symbiotic relationship between the bacteria and their host. Key enzymatic genes, RubisCO, and ATP citrate lyase (ACL) type II associated with these carbon fixation cycles, were identified as proteins in Lamelibrachia symbionts (Table S#). Our results further corroborate the idea that the presence of rTCA in addition to CBB pathways for carbon fixation might be common in all vent-dwelling and seep-living vestimentiferan endosymbionts. We also identified several key components related sulfide and nitrogen metabolic pathways, and the results are largely consistent with previous analyses.

Hemoglobin evolution

Figure 3. Hemoglobin gene family in Lamellibrachai genome. (A) Partial alignment of sampled siboglinid HB subunit A1, A2, B1, B2 sequences. Only positions contains free cysteines or cysteine residues are involved in disulfid bridge were shown in the alignment. (B) Gene tree of siboglinid HB subunit A1, A2, B1 and B2. Only Lamellibrachia* HB sequences (from NCBI or this study) were labeled at the tips. (C) Lamellibrachia HB sequences identifed that are highly expressed in the trophosome tissue or as proteins from proteomic data.

Partial alignment and gene tree of siboglinid HBs is shown in Figure #. The expression of Hb subunit A1, A2, B1, B2, and linker L1-L4 were among the highest in the trophosome compared to the plume and vestimentum, and we were able to identify most of them in the mass spectrometry. Consistent with previous studies (Zal et al. 1996), a single copy of A2 and B2 HB was identified in all siboglinids with a conserved-free cystine (i.e., cysteine residues not involved in disulfide bridges) at position 77 and 67, respectively. Except for A2 and B2 HB globins in Lumbricus terrestris, Cystine residues in same positions are also identified in Cirratulus spectabilis, Sabella pacifica , and Sternapsis sp. from sulfide-free and Arenicola marina living in sulfide-rich environments (Fig. 3). These results support the hypothesis that free cysteine residues in A2 and B2 globins have been conserved before the radiation of the Siboglinidae and potentially involved in H2S detoxification process (Bailly et al. 2002).

Surprisingly, we found a significant expansion of B1 globins in Lamellibrachia genome with 25 copies. Riftia might also contain more copies since three B1 globins were identified using a set of degenerate primers by PCR in the previous study (Bailly et al. 2002). Noticeably, we found that eight copies of “Lamellibrachia” B1 sequences also contains a free cysteine at position 77, same position as free cysteine in A2 globins. Five out of eight copies were highly expressed in the trophosome, and one copy was identified at the protein level. Unlike A2 and B2 globins, it was long thought that B1 globins can only bind O2 and lack of the ability of sulfid binding. Although the presence of free cysteine residues is not proof of sulfide binding also occurs in B1 globins, it is at least possible that the hemoglobin structure and sulfide binding function is more complex than previous thought due to unexpected large expansion and high expression of “Lamellibrachia” B1 sequences.

Instead of free cysteines, another hypothesis suggested that the H2S binding affinity has been mediated by the zinc moieties bound to amino acid residues at the interface between pairs of A2 chains in Riftia (Flores et al. 2005). Figure S# shows an alignment of amino acid sequences of A2 globins from Siboglinidae and other annelids. Amino acids involved in the binding to of zinc ions in Riftia A2 globin pairs are highlighted. Each A2 chain contains one complete Zn2+-binding site and shares two others with adjoining A2 subunits. The Zn2+-binding site contained within A2 chain is composed of three His residues (B12, B16, and G9) (Flores et al. 2005). However, none of these sites are conserved across siboglinids, even in vestimentiferans (Fig. S#). For example, B16 His site is replaced by a variety of amino acids in seep-living vestimentiferans. Thus, it is less likely that the zinc sulfide binding mechanism is a conserved function in siboglinids, and might only be a derived function just for vent-living vestimentiferans.

Immunity function

Figure 4. Putative Toll-like recepor pathway in Lamellibrachia genome. (A) Putative TLR4-MYD88 and NF-kB pathway in Lamellibrachia might be essential for immunity and response to symbionts and pathogens. (B) Toll-like recepor gene trees from selected lophotrochozoan genomes.

The immune function between hosts and symbionts is a key evolutionary driver that is under strong selective pressure and has potential implications in aging. However, the genetic machinery and functionality of the immune system in most chemosynthetic symbioses have not been extensively characterized. Our detailed analysis revealed several genes related to immune function found duplications in Lamellibrachia compared with other lophotrochozoans (Table S#). Toll-like receptor (TLRs) provides vital cellular and molecular defense against invading pathogens and recognition of host-microbial symbiosis (Chu and Mazmanian 2013). Consistent with previous analysis (Luo et al. 2018), we found that TLR genes were generally expanded in most lophotrochozoans, especially in bivalves (Fig. S#). Noticeably, we found several lineage-specific expansion of putative TLR genes in Lamellibrachia genome with 31 copies, whereas only five copies were found in Capitella genome, suggesting TLR genes in Lamellibrachia may play a more relevant role for symbiosis.

The initial physical encounter between tubeworms and symbionts occurs in an extracellular mucous secreted by pyriform glands by newly settled larvae (Nussbaumer, Fisher, and Bright 2006). In mucus matrices, symbionts can attach to the host using extracellular components secreted from symbionts, such as lipopolysaccharide - a major component of the outer membrane of Riftia symbionts. TLR4 orthologs have been hypothesized to recognize lipopolysaccharide in several studies. Although our TLR tree was not fully resolved for TLR orthology assignment based on humans TLR orthologs potentially due to lack of ecdysozoan and most deuterostomic taxa, we found a genomic expansion of TLR4-like genes in Lamellibrachia (7 copies) compared with Capitella (1 copy) from KEGG mapping. Moreover, the infection process finishes simultaneously with massive apoptosis of skin tissue during symbionts travel from host epidermal cells into trophosome (Nussbaumer, Fisher, and Bright 2006). This result is consistent with that recognition of lipopolysaccharide (LPS) by TLR4 can result in induction of signaling cascades that lead to activation of NF-kb and the production of proinflammatory cytokines (Chu and Mazmanian 2013). Interestingly, we also found genomic expansions of Tumor necrosis factor receptor (TRFR) and TNF recepotr associated factor (TNAF) (Table S#), both play vital roles in the NF-kb pathway. Although it is still not clear how the hosts distinguish between symbionts and pathogens in most symbioses, alkaline phosphatase has been shown to be involved in the maintenance of hpmeostasis of interesting commensal bacteria in the squids, mouse, and zebrafish. The commensal bactieria-derived LPS signaling via TLR4-MYD88 leads to upregulation of intestinal alkaline phosphatase and prevents inflammatory responses to the resident microbiota. Importantly, we also identified eight copies of alkaline phosphatase in Lamellibrachia, whereas only one copy was found in Capitella, suggesting a potential mechanism of in response to symbionts and pathogens, thus facilitating colonization by the endosymbionts. Although further analysis is warranted to investigate these hypotheses explicitly, TLR4-like/MYD88 and NF-kB pathway might be essential to host immunity and in response to symbionts and pathogens.

Besides TLRs, we also found evidence of the genomic expansion of NLRP family, which play a key role in the inflammasome, an innate immune system protein complex that recognizes infectious pathogens and regulates the activation of inflammatory caspases. Moreover, We found large expansions of several members of Sushi domain proteins that could be involved recognition and adhesion between hosts and symbionts.

Aging and cell cycle

The extreme longevity of seep-living vestimentiferans is expected to involve genes affecting aging. However, genetic mechanisms involved in aging on Lamellibrachia are incredibly scarce. Besides genomic expansions of genes related to innate immunity mentioned above, our analysis highlighted several genes in Lamellibrachia that may play a direct role in aging (Table #). We found evidence of positive selection of RNA polymerases, includes RNA polymerase I (RNAPI), RNA polymerase II (RNAPII) and RNA polymerase III subunit K (RNAPIIIK). Importantly, RNAPII can be used for recognition proteins in cellular responses to DNA damage. RNAPII can stall at a damaged DNA base may constitute the most specific signal for DNA repair, DNA damage checkpoints and apoptosis (Lindsey-Boltz and Sancar 2007) (Derheimer et al. 2007). The previous study has also suggested that RNAPII and RNAPIII might involve in DNA repair for UV damage in yeasts(Charton et al. 2015).

Moreover, we identified the signature of positive selection in several Histone H2 members and Lysine-specific histone demethylase 1 (Table #). Changes of histones can lead to profound chromatin structure alterations thereby influencing a multitude of biological processes ranging from transcriptional regulation to genome stability (Bönisch and Hake 2012). Changes of specific histone methyltransferases and demethylases and alterations in histone methylation could have a huge impact on lifespan and aging phenotypes, as seen in yeast, Caenorhabditis elegans and mammals (see review in (McCauley and Dang 2014)). These results suggest that RNA polymerases regulate the strength of DNA repair mechanisms and histone modifications or methylations may play a critical role in successful aging in Lamellibrachia.

We also found selection signals in microtubulin structures (tubulin alpha, beta and gamma chains), suggesting that changes in cytoskeletal dynamics may have been important in the evolution of Lamellibrachia. It has been shown that cell proliferation activities of the seep-living Lamellibrachia and vent-living Rifita are higher than in any other characterized invertebrates, being only comparable with tumor and wound-healing processes (Pflugfelder, Cary, and Bright 2009). Microtubules play an essential role in cell division, and expression of tubulin content increases as the cells progress toward mitosis (???). Interestingly, tubulins were among the strongest selection genes in long-lived vertebrate giant tortoise genome (Quesada et al. 2018), suggesting that alterations of tubulins might be related to the increase in the number of cellular divisions associated to gigantism and high proliferation rates in Lamellibrachia.

Phylogeny and molecular clock of tubeworm

Figure 5. Molecular clock tree of tubeworm origin.

Two additional Osedax taxa were added in the study compared to the previous siboglinid phylogenomic analyses (Li et al. 2017) (Table #). The final supermatrix dataset contains 191 genes single-copy orthologs. Bayesian inference with a relaxed molecular clock (Fig. 2) recovered the same topology as ML analysis in IQTree with strong nodal support (Fig. #). Both analyses strongly support Osedax is most closely related to the Vestimentifera + Sclerolinum clade and Frenulata is the early diverging group, as recently reported (Li et al. 2015). Within Vestimentifera, Lamellibrachia is sister to the remaining sampled vestimentiferans. Molecular clock analyses based on phylogenomic dataset suggest modern siboglinid diversity originated in Mesozoic (223MYA ± 80 MY), conflicting with previous hypotheses indicate a Late Mesozoic or Cenozoic approximately 50-126 MYA. The previous analyses are solely based on the estimation of nucleotide divergence of COI sequences on limited taxa sampling (mainly vestimentiferans). However, mitochondrial genes of vestimentiferans may have experienced a “slow-down” in rate of nucleotide subsitution relative to other siboglinid lineages (Li et al. 2015).

Moreover, molecular clock analysis suggests a younger split of Vestimentifernas during the Cenozoic (60 MYA). Recent analyses indicated that Jurassic tubes from Figueroa deposits are likely to have been made by vestimentiferans (Georgieva et al. 2017), but this hypothesisis is not supported in our analysis. Thus, we provide the molecular clock of the siboglinid phylogenetic tree, placing a common siboglinid ancestor as far back as the Jurassic. Siboglinid vestimentiferans might exploit and adapt to vents and seeps during Cenozoic (<65 MYA). This adds to the growing evidence that the Cenozoic was a key period for the radiation of most dominant invertebrate taxa now occupying in deep-sea chemosynthetic communities (Vrijenhoek 2013).